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Sir: 

Prior to examination of the above-referenced application, please enter the 
following amendments and remarks. 

IN THE SPECIFICATION: 



Please replace the paragraph beginning on page 2, line 14, with the 
following rewritten paragraph: 



The field of molecular modeling has successfully simulated the motion 
(molecular dynamics or MD) and determined energy minima or rest states (static 
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analysis) of many complex molecular systems by computers. Typical molecular 
modeling applications have included enzyme-ligand docking, molecular diffusion, 
reaction pathways, phase transitions, and protein folding studies. Researchers in the 
biological sciences and the pharmaceutical, polymer, and chemical industries are 
beginning to use these techniques to understand the nature of chemical processes in 
complex molecules and to design new drugs and materials accordingly. Naturally, the 
acceptance of these tools is based on several factors, including the accuracy of the results 
in representing reality, the size and complexity of the molecular systems that can be 
modeled, and the speed by which the solutions are obtained. Accuracy of many 
computations has been compared to experiment and generally found to be adequate 
within specified bounds. However, the use of these tools in the prior art has required 
enormous computing power to model molecules or molecular systems of even modest 
size to obtain molecular time histories of sufficient length to be useful. 

Please replace the paragraph beginning on page 3, line 6, with the 
following rewritten paragraph: 



Substantial work has been completed in reducing the computational load 
for molecular models, such as the reduction of model complexity by constraining higher 
order modes with rigid body assumptions, simplifying the model with rigid or flexible 
substructuring, Order(N) dynamics, efficient implicit solvent models, and multipole 
methods for the force field models (see, for example, U.S. Patent No. 5,424,963 on the 
commercial MBO(N)D software package). Co-pending applications, US. Appln. No. 
10/053,253, entitled "METHOD FOR LARGE TMESTEPS IN MOLECULAR 
MODELING," and U.S. Appln. No. 10/053,354, entitled "METHOD FOR RESIDUAL 
FORM IN MOLECULAR MODELING" both of which are filed of even date, claim 
priority from the previously cited provisional patent applications and which are assigned 
to the present assignee, and which are incorporated by reference herein, describe further 
improvements in molecular models and numerical methods. 
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Please replace the paragraph beginning on page 3, line 19, with the 
following rewritten paragraph: 

The primary reason molecular simulations are so slow is that currently 
used or applied numerical methods require very small timesteps, typically between 1 and 
10 femtoseconds (10 -15 to 10" 14 seconds). Each timestep taken requires the computation 
of a new state (position and motion for each atom) of the particular molecular model, and 
then computation of the new set of forces resulting from the new state. For example, 
molecular dynamics simulations of the complex behavior of large molecules, such as the 
folding of a protein, typically need to cover a time span from at least a microsecond up to 
several seconds or even minutes. With techniques currently in common use, this results 
in the requirement to take 10 9 to 10 16 timesteps in the computer simulation. The per-step 
computations of the state, and especially the forces, grow very expensive as the problem 
size increases. Even with the fastest computers available today, months, years or even 
centuries of computer time are required to solve such problems even for systems of 
modest size. 

Please replace the paragraph beginning on page 4, line 6, with the 
following rewritten paragraph: 

This common-sense belief is incorrect, however. The computer science 
sub-discipline of numerical analysis has produced an extensive theory of numerical 
integration for problems in which high frequencies exist but are of little interest. These 
problems are termed "stiff problems (see, for example, Hairer and Wanner, Solving 
Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2 nd ed. , 
Springer, 1996). In these cases, it is the stability of the integration method, not the 
required solution accuracy, which limits the timestep. Integration methods exist which 
have unconditional stability, meaning that in theory they can take arbitrarily large 
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timesteps. Such methods have a mathematical property called "L-stability." Only so- 
called "implicit" integration methods can be L-stable, but very few implicit integration 
methods actually are L-stable. Previously cited co-pending U.S. Patent Appln, No. 
10/053,253, entitled "METHOD FOR LARGE TIMESTEPS IN MOLECULAR 
MODELING," covers the use of implicit and in particular L-stable integration methods. 



Please replace the paragraph beginning on page 4, line 24, with the 
following rewritten paragraph: 

The same problem of high frequency molecular vibration for MD 
simulations causes problems for the calculation of Jacobians. For example, the repulsive 
van der Waal's forces that are generated as the electron orbitals of two atoms begin to 
overlap must be accounted for in a molecule. This quantum mechanical effect is often 
approximated by the so-called Lennard Jones potential (Rapaport, op. cit.) 9 which models 
the repulsive force as being proportional to 1 / r 13 , where r is the distance between the 
centers of the atoms. This is an extremely stiff interatomic force which is characteristic 
of molecular dynamics (MD) simulations and poses particular difficulty for any 
numerical integration scheme used to advance time during the simulation. As a result 
and as mentioned previously, most prior art MD integration schemes have timesteps 
limited by the high frequencies generated by these extremely stiff repulsive forces. And 
in particular, the stiffness of the atomic forces greatly magnify the difficulty of forming 
certain required Jacobian matrices. Such Jacobian matrices are a necessary ingredient of 
any stable implicit integration scheme, such as described in the immediately cited co- 
pending application. 

Please replace the paragraph beginning on page 5, line 24, with the 
following rewritten paragraph: 
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All general-purpose simulation codes provide routines to compute 
Jacobians numerically as follows. For a given equation to integrate y = f (y, t) , the 
desired Jacobian is J = df Idy and is numerically computed: 

Ay 

where 

A/ = /U 2 -/ 
Ay = y 2 - yi 

Note that the perturbation Ay has to be selected to give a good result and 

may be difficult to choose, especially for stiff systems. In contrast, analytic Jacobians are 
computed by solving directly, or in this case algorithmically, for the equation of the 
desired derivatives. 



Please replace the paragraph beginning on page 6, line 8 , with the 
following rewritten paragraph: 

The present invention provides for a method of modeling the behavior of a 
molecule. The method has the steps of selecting a torsion angle, rigid multibody model 
for the molecule, the model having equations of motion; selecting an implicit integrator; 
and generating an analytic Jacobian for the implicit integrator to integrate the equations 
of motion so as to obtain calculations of the behavior of the molecule. The analytic 
Jacobian J is derived from an analytic Jacobian of the Residual Form of the equations of 
motion and is described as: 





dq^ 






dq 


du 




'■>„ 


du 


du 






K dq 


du j 







Juu J 



; and 
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qq dq dq qu du 

uq dq dq uu du du 

where q are the generalized coordinates, u are the generalized speeds, Wis 
a joint map matrix and M is the mass matrix and p u is the dynamic residual of the 
equations of motion, and z is -M~ l p u w, 0) . The method can also be used for any 
physical system which can be modeled by a torsion angle, rigid multibody system. 



Please replace the paragraph beginning on page 9, line 1 , with the 
following rewritten paragraph: 

The description of the preferred embodiment is divided into several 
sections. The first set of sections describes the MD simulation architecture, the 
multibody system (MBS) definitions, and Residual Form of the MBS equations for the 
subsequent descriptions. The next set of sections discusses the definition of the Jacobian, 
its role in the implicit integration method, and the computation of the analytic Jacobian 
using the Residual Form. Also shown is the superior accuracy and performance of the 
analytic Jacobian vs. the numerical Jacobian. Further efficiencies in the computation of 
the analytic Jacobian are discussed, specifically, exploiting the rigid body MBS to 
"contract" the size of the Jacobian matrix, and exploiting the topological structure of the 
MBS to eliminate unnecessary computations. 

Please replace the paragraph beginning on page 9, line 29 , with the 
following rewritten paragraph: 

In accordance with the present invention, a method for the solution of the 
equations of a molecular system is expressed in Residual Form to bypass the customary 
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step of producing the state derivatives directly. The Residual Form method has the 
following steps: 

1) Discretization of the solution variables. The specific form of discretization is 
dictated by the particular implicit integration method used to advance the 
molecular model in time. Implicit integration follows from the Residual Form. 
Implicit integration, especially L-stable integrators and other highly stable 
integrators, such as implicit Euler, RadauS, SDIRK3, SDIRK4, other implicit 
Runge-Kutta methods, and DASSL or other implicit multistep methods, also 
provide other advantages for molecular modeling. See, for example, the above- 
cited U.S. Patent Appln. No. 10/053,253, entitled "METHOD FOR LARGE 
TIMESTEPS IN MOLECULAR MODELING," filed of even date. As a 
particularly simple example, when used with implicit Euler integration, the 
discretization is as follows: 

where h is the timestep. 

Please replace the paragraph beginning on page 11, line 15, with the 
following rewritten paragraph: 

The general system architecture 48 of the software and some of its 
processes for modeling molecules in accordance with the present invention are illustrated 
in Fig. 1 . Each large rectangular block represents a software module and arrows 
represents information which passes between the software modules. The software system 
architecture has a modeler module 50, a biochemistry components module 52, a physical 
model module 54, an analysis module 56 and a visualization module 58. The details of 
some of these modules are described below; other modules are available to the public. 
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Please replace the paragraph beginning on page 1 1, line 23, with the 
following rewritten paragraph: 



The modeler module 50 provides an interface for the user to enter the 
physical parameters which define a particular molecular system. The interface may have 
a graphical or data file input (or both). The biochemistry components module 52 
translates the modeler input for a particular mathematical model of the molecular system 
and is divided into translation submodules 60, 62 and 64 for mathematical modeling the 
molecule(s), the force fields and the solvent respectively of the system being modeled. 
There are several modeler and biochemistry components modules available including, for 
example, Tinker (Jay Ponder, TINKER User's Guide , Version 3.8, October 2000, 
Washington University, St. Louis, MO). 

Please replace the paragraph beginning on page 1 1, line 32, with the 
following rewritten paragraph: 

With the translated physical parameters from the biochemistry 
components module 52, the physical model module 54 defines the molecular system 
mathematically. At the core of the module 54 is a multibody system submodule 66. The 
analysis module 56, which communicates with the physical model module 54 and the 
visualization module 58, provides solutions to the computational models of the molecular 
systems defined by the physical model module 54. The analysis module 56 consists of a 
set of integrator submodules 68 which integrate the differential equations of the physical 
model module 54. The integrator submodules 68 advance the molecular system through 
time and also provide for static analyses used in determining the minimum energy 
configuration of the molecular system. The analysis module 56 and its integrator 
submodules 68 contain most of the subject matter of the present invention and are 
described in detail below. 
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Please replace the paragraph beginning on page 12, line 11, with the 
following rewritten paragraph: 



The visualization module 58 receives input information from the 
biochemistry components module 52 and the analysis module 56 to provide the user with 
a three-dimensional graphical representation of the molecular system and the solutions 
obtained for the molecular system. Many visualization modules are presently available, 
an example being VMD (A. Dalke, et al. 9 VMD User's Guide , Version 1.5, June 2000, 
Theoretical Biophysics Group, University of Illinois, Urbana, Illinois). 

Please replace the paragraph beginning on page 12, line 3 1 , with the 
following rewritten paragraph: 

The MBS is an abstraction of the atoms and effectively rigid bonds that 
make up the molecular system being modeled and is selected to simplify the actual 
physical system, the molecule in its environment, without losing the features important to 
the problem being addressed by the simulation. With respect to the general system 
architecture illustrated in Fig. 1, the MBS does not include the electrostatic charge or 
other energetic interactions between atoms nor the model of the solvent in which the 
molecules are immersed. The force fields are modeled in the submodule 62 and the 
solvent in the submodule 64 in the biochemistry components module 52. 

Please replace the paragraph beginning on page 13, line 31 , with the 
following rewritten paragraph: 

The symbol for the vector from one point to another contains the name of 
the two points. Thus, r PQ is the vector from the point P to point Q. A vector 
representing the velocity of a point in a reference frame contains the name of the point 
and the reference frame: N v p . Certain symbols to be introduced later relate two reference 
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frames. In this case, the symbol contains the name of two frames. Thus, l C k is the 
direction cosine matrix for the orientation of frame k in frame L This symbol refers to the 
direction cosine matrix for a typical body in its parent frame. Thus, l C k (j) indicates the 
direction cosine matrix of the actual body j in question. The left and right superscripts do 
not change with the body index. This is also true for the other symbols. An asterisk 
indicates the transpose: H*(k) , for example. A tilde over a vector indicates a 3 by 3 
skew-symmetric cross product matrix: vw = v x w . is an i by i identity matrix., and 
Qi is a zero vector of length i and 0. is an i by i zero matrix. 



Please replace the paragraph beginning on page 20, line 2 , with the 
following rewritten paragraph: 

With the first kinematic calculations described above, the residual 
computation for the Residual Form method can be determined. A detailed description of 
the Residual Form and its application to molecular modeling is found in the previously 
cited co-pending U.S. Patent Appln. No. 10/053,354, entitled, "METHOD FOR 
RESIDUAL FORM IN MOLECULAR MODELING," filed of even date. This 



computation fills in two partitions of the vector 



given previously. The first 



\PuJ 

partition is called p q , the kinematic residual, and the second partition is called p u , the 
dynamic residual. The kinematic residual is computed from the difference between a q , 
which is passed-in from the (implicit) integration submodules 66, and the derivative 
computed by each joint: 

q-W(q)u = p q 



Please replace the paragraph beginning on page 23, line 24 , with the 
following rewritten paragraph: 
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An important ingredient of this integration process is formation of the 
Jacobian of the differential equations. This is 

dy 

Since the function /is itself computed by an algorithm rather than by an explicit formula, 
the Jacobian computation represents a substantial amount of work. In the simplest 
approach, the Jacobian can be formed numerically by differencing the derivative routine. 
This is a delicate operation because the quality of the Jacobian is a tradeoff between 
round-off and truncation errors. Typically half the working precision in the result is 
obtained by choosing a good perturbation size in the difference scheme. In practice, 
though, this is difficult to do. 

Please replace the paragraph beginning on page 24, line 10 , with the 
following rewritten paragraph: 

In general, G, the iteration matrix used in the Newton loop of the implicit 
integrator has the form G-E-aJ, where E is the identity matrix and a is some 
scalar function of the timestep. See the previously referenced U.S. Patent Appln. No. 
10/053,253, entitled "METHOD FOR LARGE TIMESTEPS IN MOLECULAR 
MODELING," filed of even date, for a description of implicit integrators. Changes in 
step size require refactoring G , but not reforming J . Reforming J is needed only when 
the Jacobian is needed at a new state. G is used in a linear subproblem within a Newton 
loop. The following is solved: 

where r(y) is the residual function for that particular implicit integration method. 
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Please replace the paragraph beginning on page 27, line 16 , with the 
following rewritten paragraph: 



upon two kinds offerees: 'motion forces' and external forces. The motion forces are 
computed directly by the multibody system. The external forces are available to the 
multibody system from a force modeling routine that computes the various interatomic 
forces such as electrostatics and solvents. A similar procedure is followed when 
computing the Jacobian. The multibody system builds the Jacobian of the motion forces, 
and combines it with the Jacobian of the external forces. 



their intrinsic coordinates are brought into the MBS coordinates. This requires the use of 



the chain-rule to transform between intrinsic and the MBS generalized coordinates. It is 
important that each effect co-computes its function value and Jacobian, because many of 
the same terms are needed for each computation. Each effect is transformed into a set of 
spatial loads T ejTect (k) , where k is the index of a generic body in the system. The totality 

of these effects is given the symbol T(k). 

Please replace the paragraph beginning on page 34, line 3, with the 
following rewritten paragraph: 

Note that the row-reduction procedure only needs to be done once before 
computing the Residual Jacobian. The overhead of performing the reduction is more than 
offset by the reduced cost of the smaller matrix vector products which must be formed. 



At a high level, the residual computation can be considered to depend 




following rewritten paragraph: 



Please replace the paragraph beginning on page 27, line 31, with the 



The same exchange occurs to compute Jacobians. The native Jacobians in 



i 



I W illi ! 
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Note that in forming ^ fa Qre ^ n0 nee( j t0 save ^ e e } emen t s 0 f the atomic Force 

dr(p,s) 

Jacobian. That is, each element ^Htlll only needs to be available while its contribution 

dr(p,s) 

to is being computed. So, more than one element of the big Force Jacobian is 

dr(p,s) 

not required at a time. 



Please replace the paragraph beginning on page 37, line 17, with the 
following rewritten paragraph: 

1. Given (q,u) compute z = -M~ 1 /? w (#,h,0) using the Direct Form method. Then 
set ii = z and recompute , then p u , which recomputes r(fc) . 

Please replace the paragraph beginning on page 38, line 17, with the 
following rewritten paragraph: 

Each interbody direction cosine matrix (and all the joint-specific) quantities 
depend only on the generalized coordinates of an individual joint. Thus, 
l dC k (k) is nonzero only when the derivative is taken with respect to any of the 
coordinates for body k. To properly 'seed' the recursions being studied, a vector 
dq is passed in to the routine. For Jacobian computation we set one entry is set to 

1 , and all the other entries to 0. Then, the needed preliminary quantities are 
generated by a typical loop: 

dq{k) 

The partial derivatives of the direction cosine matrices are generated analytically 
and displayed in the section, "Interface Contraction" above. These partial 
derivatives do not depend upon the particular column of the Jacobian that is being 



m 
bi 
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computed. Setting a particular entry of dq to 1, and all the rest to 0 has the effect 
of annihilating the correct subset of the seed quantities. 
dr &Qt (k) 



■ , the partial derivative of the interbody spanning vector is given by 



Br QA {k) 
dq k 

dr Q& {k) 
dq k 

8r Q& (k) 
dq k 

dr m (k) 
d 1k 



= A{k), slider 



= 1U ball 



= 0,, pin 

0 0 0 0 1 0 0" 

0 0 0 0 0 1 0 

0 0 0 0 0 0 1 



free 



Z{k) here refers to the body's sliding axis which connects it to its parent body. 
dH(k) 



-, the partial derivative of the joint map is 



dH(k) n . . ... 

=0 6, pin, slider 

dq k 

f^=[a a], ban 

dq k L_ _ J 



dH{k) 



a 



a 



d'Cjk) 



for £ = 1,2,3,4, free 



^^ = 0 6 for k = 5,6,7, free 



Please replace the paragraph beginning on page 39, line 13, with the 
following rewritten paragraph: 
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Since arbitrary perturbations to a set of Euler parameters do not produce a 
pure rotation, column contraction cannot be used when computing the 
corresponding column of the Jacobian. The row-reduced Force Jacobian can still 
(and must) be used. 

Please replace the paragraph beginning on page 42, line 13 , with the 
following rewritten paragraph: 



10. Back-solve the result of previous step with the mass-matrix to 
dq 

obtain the desired — : 
dq 

dq dq 

The back-solve operation is accomplished in the Direct Form method routine by 
processing a residual vector into a u vector. The Second Kinematics Calculations only 
needs to be performed once for the whole Jacobian, since the back-solves are done at the 
nominal value of the state. In fact, the Second Kinematics routine must have been called 
in Step 1 while computing z, so the variables should still be cached. 

Please replace the paragraph beginning on page 44, line 7, with the 
following rewritten paragraph: 

14. Back-solve the result of previous step with the mass-matrix to 
du 

obtain the desired — : 
du 

du du 

The back-solve operation is accomplished in the Direct Method routine by processing a 
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residual vector into a u vector. The Second Kinematics Calculations only need to be 
performed once, since the back-solves are done at the nominal value of the state. In fact, 
the Second Kinematics routine must have been called in Step 1 while computing z, so the 
variables should still be cached. 



Please replace the paragraph beginning on page 47, line 2 1 , with the 
following rewritten paragraph: 

Therefore, while the foregoing is a complete description of some of the 
embodiments of the invention, it should be evident that various modifications, 
alternatives and equivalents may be made and used. Accordingly, the above description 
should not be taken as limiting the scope of the invention which is defined by the metes 
and bounds of the appended claims. 

IN THE CLAIMS: 



Please add new claims 17-44 as follows: 



1 7. A method of modeling the behavior of a molecule, comprising 
selecting a torsion angle, rigid multibody model for said molecule; and 
generating a partition of an analytic Jacobian J for said model, said 
partition selected from the group comprising J qq9 J qu , J m and J uu , and 



8q 


dq^ 






dq 


du 




(J J \ 

J qq J qu 


du 


du 
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du j 
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where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and it are derivatives of said 
generalized speeds with respect to time. 

18. The method of claim 1 7 wherein said model has equations of 

motion and 

qq dq dq qu du 

j =* = ^-iSflto) and j m = *=-M-*&t^ 
dq dq m du du 

where Wis a joint map matrix defining relationships between q and u, M 
is the system mass matrix defining relationships between u and the generalized forces on 
said bodies, p u is the dynamic residual of said equations of motion, and z is 
~M> a ( 9 ,w,0). 

1 9. The method of claim 1 8 wherein said molecule comprises a 
number of atoms and wherein said selecting step further comprises 

reducing said number of atoms to a smaller number of rigid bodies so as to 
reduce the number of computations for Jacobian partitions J uq and J uu . 

20. The method of claim 19 wherein said number of computations is 
reduced approximately by the average number of atoms in each of said rigid bodies. 



comprising 



21. A method of modeling the behavior of a physical system, 
selecting a torsion angle, rigid multibody model for said physical; and 
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generating a partition of an analytic Jacobian J for said physical model, 
said partition selected from the group comprising J qq , J qu , J uq and J uu , and 



(dq dq} 

dq du 

du du 

{dq du) 



where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and u are derivatives of said 
generalized speeds with respect to time. 



22. The method of claim 2 1 wherein said physical model has equations 



of motion and 



= dq = d{Wu) 
qq dq dq 



and J =^- = W 
qu du 



where Wis a joint map matrix defining relationships between q and w, M 
is the system mass matrix defining relationships between u and the generalized forces on 
said bodies, p u is the dynamic residual of said equations of motion, and z is 



23. Computer code for simulating the behavior of a molecule, said 
code comprising 

a first module for a torsion angle, rigid multibody model for said 
molecule; and 
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a second module for generating a partition of an analytic Jacobian J for 
said model, said partition selected from the group comprising J qq , J qu , J m and«/„„, and 



dq 


dq 


dq 


du 


du 


du 


dq 


du 



(J J 

J qq J qu 



\*^uq 



J. 



where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and u are derivatives of said 
generalized speeds with respect to time. 



24. The computer code of claim 23 wherein said model has equations 



of motion, and 



= dq = d{Wu) 
qq " dq dq 



09 dq dq m du du 

where Wis a joint map matrix defining relationships between q and «, M 
is the system mass matrix defining relationships between u and the generalized forces on 
said bodies, p u is the dynamic residual of said equations of motion, and z is 
-M> M (?,«,<)). 



25. The computer code of claim 24 wherein said molecule comprises a 
number of atoms and wherein said first module further comprises 

code reducing said number of atoms to a smaller number of rigid bodies so 
as to reduce the number of computations for Jacobian partitions J m and J m « 
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26. The computer code of claim 25 wherein said number of 
computations is reduced approximately by the average number of atoms in each of said 
rigid bodies. 

27. Computer code for simulating the behavior of a physical system, 
said code comprising 

a first module for a torsion angle, rigid multibody model for said physical 

system; and 

a second module for generating a partition of an analytic Jacobian J fox 
said model, said partition selected from the group comprising J qq9 J qu , J uq and J uu , and 



J = 



'dq dq} 

dq du 

du du 

K dq dU; 



J. 



uu ) 



where q are generalized coordinates, q are derivatives of said generalized 
coordinates with respect to time, u are generalized speeds, and u are derivatives of said 
generalized speeds with respect to time. 



28. The computer code of claim 27 wherein said model has equations 



of motion, and 



j JA= £&) and 

" dq dq 



qu du 



" dq dq du du 



where FT is a joint map matrix defining relationships between q and «, M 
is the system mass matrix defining relationships between u and the generalized forces on 
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said bodies, p u is the dynamic residual of said equations of motion, and z is 
-M- l p u (q 9 u,0) 

29. The method of claim 1 further comprising calculating said 
analytical Jacobian in no more than order N p_l computations, where N is proportional to 
the number of rigid bodies in said model, and wherein a calculation of a numerical 
Jacobian for said model requires order N p computations. 

30. The method of claim 5 further comprising calculating said 
analytical Jacobian in no more than order N p_i computations, where N is proportional to 
the number of rigid bodies in said model, and wherein a calculation of a numerical 
Jacobian for said model requires order N p computations. 

3 1 . The computer code of claim 9 wherein said second module further 
comprises code for calculating said analytical Jacobian in no more than order N p_1 
computations, where N is proportional to the number of rigid bodies in said model; and 
wherein a calculation of a numerical Jacobian for said model requires order N p 
computations. 

32. The computer code of claim 13 wherein said second module 
further comprises code for calculating said analytical Jacobian in no more than order N p ~ l 
computations, where N is proportional to the number of rigid bodies in said model; and 
wherein a calculation of a numerical Jacobian for said model requires order N p 
computations. 

33. The method of claim 1 wherein a calculation of a numerical 
Jacobian for said model has an accuracy of K digits and further comprising 

calculating said analytical Jacobian for said model with an accuracy of at 

least 2K digits. 

34. The method of claim 5 wherein a calculation of a numerical 
Jacobian for said model has an accuracy of K digits and further comprising 
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calculating said analytical Jacobian for said model with an accuracy of at 

least 2K digits. 

35. The computer code of claim 9 wherein a calculation of a numerical 
Jacobian for said model has an accuracy of K digits and said second module further 
comprises 

code for calculating said analytical Jacobian for said model with an 
accuracy of at least 2K digits. 



numerical Jacobian for said model has an accuracy of K digits and said second module 
further comprises 

code for calculating said analytical Jacobian for said model with an 
accuracy of at least 2K digits. 

37. The method of claim 1 wherein said analytic Jacobian is derived 
from an analytic Jacobian of the Direct Form of the equations of motion. 

38. The method of claim 37 wherein said analytic Jacobian / 

comprises 



36. The computer code of claim 13 wherein a calculation of a 
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where q are the generalized coordinates, u are the generalized speeds, Wis 
a joint map matrix defining relationships between q and u, M~ l is the operational inverse 
of the system mass matrix M which defines relationships between u and the generalized 
forces/ on said bodies. 

39. The method of claim 5 wherein said analytic Jacobian is derived 
from an analytic Jacobian of the Direct Form of the equations of motion. 



40. The method of claim 39 wherein said analytic Jacobian J 



comprises 



' qq 



( dq 


dq} 




dq 


du 




du 


du 




K dq 


duj 




dq_ 


d(Wu) 


dq 


dq 



(J., J, 



qu 



; and 



qu 



du 



■ W 



dq m du 



where q are the generalized coordinates, u are the generalized speeds, W is 
a joint map matrix defining relationships between q and w, M~ l is the operational inverse 

of the system mass matrix M which defines relationships between u and the generalized 
forces / on said bodies. 

4 1 . The computer code of claim 9 wherein said analytic Jacobian is 
derived from an analytic Jacobian of the Direct Form of the equations of motion. 



42. The computer code of claim 41 wherein said analytic Jacobian J 

comprises 



Dan E. Rosenthal 
Application No.: 10/053,348 
Page 24 



PATENT 



J = 



' dq_ df" 

dq du 

du du 

K dq du j 



(J qq J qu 



and 



= dq = d(Wu) 
qq ~ dq~ dq 



and J = — = W 

qu du 



J. 



a(M-(/)) 



uq 



dq 



and J = 



d(M~ l (f)) 



du 



where 9 are the generalized coordinates, u are the generalized speeds, J^is 
a joint map matrix defining relationships between q and u 9 M~ l is the operational inverse 
of the system mass matrix M which defines relationships between u and the generalized 
forces / on said bodies. 

43. The computer code of claim 13 wherein said analytic Jacobian is 
derived from an analytic Jacobian of the Direct Form of the equations of motion. 



44. The computer code of claim 43 wherein said analytic Jacobian J 



comprises 



f dq dq} 

dq du 

du du 

^dq du j 



; and 



j „M.»fl5l and J = *~W 
« dq dq * du 



m dq du 
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where q are the generalized coordinates, u are the generalized speeds, Wis 
a joint map matrix defining relationships between q and u, M" 1 is the operational inverse 
of the system mass matrix M which defines relationships between u and the generalized 
forces / on said bodies. 



IN THE DRAWINGS: 

Four sheets of replacement drawings are being submitted. 



REMARKS 

Consideration of the patent application, as preliminarily amended, is 
respectfully requested. 

By this amendment of the specification, the applicant has corrected 
various typographical errors, filled in the blank spaces in the text of the specification with 
now issued U.S. Application Nos., and described his invention with greater precision. 
For example, the change in the mathematical expression on page 5 reflects the better 

statement that the numerically computed — is approximately equal to the analytical 

Ay 

expression of the Jacobian J. Likewise, the changes to the equation on page 39 make that 
equation more explicit with respect to the index k. 

With respect to the claims, the applicant has also added new claims 17-44 
to more particularly point out and distinctly claim the subject matter which he regards as 
his invention entitled. Claims 1-44 are now pending. 

If the Examiner believes a telephone conference would expedite 
prosecution of this application, please telephone the undersigned at 650-326-2400. 
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VERSION WITH MARKINGS TO SHOW CHANGES MADE 

The paragraph beginning on page 2, line 14, has been amended as follows: 

The field of molecular modeling has successfully simulated the motion 
(molecular dynamics or [(] MD [)]) and determined energy minima or rest states (static 
analysis) of many complex molecular systems by computers. Typical molecular 
modeling applications have included enzyme-ligand docking, molecular diffusion, 
reaction pathways, phase transitions, and protein folding studies. Researchers in the 
biological sciences and the pharmaceutical, polymer, and chemical industries are 
beginning to use these techniques to understand the nature of chemical processes in 
complex molecules and to design new drugs and materials accordingly. Naturally, the 
acceptance of these tools is based on several factors, including the accuracy of the results 
in representing reality, the size and complexity of the molecular systems that can be 
modeled, and the speed by which the solutions are obtained. Accuracy of many 
computations has been compared to experiment and generally found to be adequate 
within specified bounds. However, the use of these tools in the prior art has required 
enormous computing power to model molecules or molecular systems of even modest 
size to obtain molecular time histories of sufficient length to be useful. 

The paragraph beginning on page 3, line 6, has been amended as follows: 

Substantial work has been completed in reducing the computational load 
for molecular models, such as the reduction of model complexity by constraining higher 
order modes with rigid body assumptions, simplifying the model with rigid or flexible 
substructuring, Order(iV) dynamics, efficient implicit solvent models, and multipole 
methods for the force field models (see, for example, U.S. Patent No. 5,424,963 on the 
commercial MBO(N)D software package). Co-pending applications, U.S. Appln. No. 
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10/053,253 . entitled "METHOD FOR LARGE TIMESTEPS IN MOLECULAR 
MODELING," and U.S. Appln. No. 10/053354 . entitled "METHOD FOR RESIDUAL 
FORM IN MOLECULAR MODELING," both of which are filed of even date, claim 
priority from the previously cited provisional patent applications and which are assigned 
to the present assignee, and which are incorporated by reference herein, describe further 
improvements in molecular models and numerical methods. 

The paragraph beginning on page 3, line 19, has been amended as follows: 

The primary reason molecular simulations are so slow is that currently 
used or applied numerical methods require very small timesteps, typically between 1 and 
10 femtoseconds (10~ i5 to 10" 14 seconds). Each timestep taken requires the computation 
of a new state (position and motion for each atom) of the particular molecular model, and 
then computation of the new set of forces resulting from the new state. For example, 
molecular dynamics simulations of the complex behavior of large molecules, such as the 
folding of a protein, typically need to cover a time span from at least a microsecond up to 
several seconds or even minutes. With techniques currently in common use, this results 
in the requirement to take 10 9 to 10 16 timesteps in the computer simulation. The per-step 
computations of the state, and especially the forces, grow very expensive as the problem 
size increases. Even with the fastest computers available today, months, years or even 
centuries of computer time are required to solve such problems even for systems of 
modest size. 

The paragraph beginning on page 4, line 6, has been amended as follows: 

This common-sense belief is incorrect, however. The computer science 
sub-discipline of numerical analysis has produced an extensive theory of numerical 
integration for problems in which high frequencies exist but are of little interest. These 
problems are termed "stiff' problems (see, for example, Hairer and Wanner, Solving 
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Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems, 2 nd ed., 
Springer, 1996). In these cases, it is the stability of the integration method, not the 
required solution accuracy, which limits the timestep. Integration methods exist which 
have unconditional stability, meaning that in theory they can take arbitrarily large 
timesteps. Such methods have a mathematical property called "L-stability ." Only so- 
called "implicit" integration methods can be L-stable, but very few implicit integration 
methods actually are L-stable. Previously cited co-pending U.S. Patent Appln. No. 
10/053.253 . entitled "METHOD FOR LARGE TIMESTEPS IN MOLECULAR 
MODELING," covers the use of implicit and in particular L-stable integration methods. 

The paragraph beginning on page 4, line 24, has been amended as follows: 

[But the] The same problem of high frequency molecular vibration for 
MD simulations causes problems for the calculation of Jacobians, For example, the 
repulsive van der Waal's forces that are generated as the electron orbitals of two atoms 
begin to overlap must be accounted for in a molecule. This quantum mechanical effect is 
often approximated by the so-called Lennard- Jones potential (Rapaport, op. cit.\ which 
models the repulsive force as being proportional to 1/r 13 , where r is the distance 
between the centers of the atoms. This is an extremely stiff interatomic force which is 
characteristic of molecular dynamics (MD) simulations and poses particular difficulty for 
any numerical integration scheme used to advance time during the simulation. As a 
result and as mentioned previously, most prior art MD integration schemes have 
timesteps limited by the high frequencies generated by these extremely stiff repulsive 
forces. And in particular, the stiffness of the atomic forces greatly magnify the difficulty 
of forming certain required Jacobian matrices. Such Jacobian matrices are a necessary 
ingredient of any stable implicit integration scheme, such as described in the immediately 
cited co-pending application. 



The paragraph beginning on page 5, line 24, has been amended as follows: 
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All general-purpose simulation codes provide routines to compute 
Jacobians numerically as follows. For a given equation to integrate y = f(y 9 t), the 

desired Jacobian is J = df I dy and is numerically computed: 



where 



1 J_ Ay 



Note that the perturbation Ay has to be selected to give a good result and 
may be difficult to choose, especially for stiff systems. In contrast, analytic Jacobians are 
computed by solving directly, or in this case algorithmically, for the equation of the 
desired derivatives. 



The paragraph beginning on page 6, line 8 , has been amended as follows: 

The present invention provides for a method of modeling the behavior of a 
molecule. The method has the steps of selecting a torsion angle, rigid multibody model 
for [said] the molecule, the model having equations of motion; selecting an implicit 
integrator; and generating an analytic Jacobian for the implicit integrator to integrate the 
equations of motion so as to obtain calculations of the behavior of the molecule. The 
analytic Jacobian J is derived from an analytic Jacobian of the Residual Form of the 
equations of motion and is described as: 



dq du 
du du 
{dq duj 



qu 



; and 
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m ~dq~ dq ""'du' du 

where q are the generalized coordinates, u are the generalized speeds, W is 
a joint map matrix and M is the mass matrix and p u is the dynamic residual of the 
equations of motion, and z is -M~ l p u (q,u, 0) . The method can also be used for any 
physical system which can be modeled by a torsion angle, rigid multibody system. 



The paragraph beginning on page 9, line 1, has been amended as follows: 



The description of the p referred embodiment is divided into several 
sections. The first set of sections describes the MD simulation architecture, the 
multibody system (MBS) definitions, and Residual Form of the MBS equations for the 
subsequent descriptions. The next set of sections discusses the definition of the Jacobian, 
its role in the implicit integration method, and the computation of the analytic Jacobian 
using the Residual Form. Also shown is the superior accuracy and performance of the 
analytic Jacobian vs. the numerical Jacobian. Further efficiencies in the computation of 
the analytic Jacobian are discussed, specifically, exploiting the rigid body MBS to 
"contract" the size of the Jacobian matrix, and exploiting the topological structure of the 
MBS to eliminate unnecessary computations. 



The paragraph beginning on page 9, line 29, has been amended as follows: 



In accordance with the present invention, a method for the solution of the 
equations of a molecular system is expressed in Residual Form to bypass the customary 
step of producing the state derivatives directly. The Residual Form method has the 
following steps: 
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2) Discretization of the solution variables. The specific form of discretization is 
dictated by the particular implicit integration method used to advance the 
molecular model in time. Implicit integration follows from the Residual Form. 
Implicit integration, especially L-stable integrators and other highly stable 
integrators, such as implicit Euler, RadauS, SDIRK3, SDIRK4, other implicit 
Runge-Kutta methods, and DASSL or other implicit multistep methods, also 
provide other advantages for molecular modeling. See, for example, the above- 
cited U.S. Patent Appln. No. 10/053.253 . entitled "METHOD FOR LARGE 
TIMESTEPS IN MOLECULAR MODELING," filed of even date. As a 
particularly simple example, when used with implicit Euler integration, the 
discretization is as follows: 

where h is the timestep. 

The paragraph beginning on page 11, line 15, has been amended as 

follows: 

The general system architecture 48 of the software and some of its 
processes for modeling molecules in accordance with the present invention are illustrated 
in Fig. 1 . Each large rectangular block represents a software module and arrows 
represents information which passes between the software modules. The software system 
architecture has a modeler module 50, a bioche mistry components module 52, a physical 
model module 54, an analysis module 56 and a visualization module 58. The details of 
some of these modules are described below; other modules are available to the public. 



The paragraph beginning on page 1 1, line 23, has been amended as 

follows: 
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The modeler module 50 provides an interface for the user to enter the 
physical parameters which define a particular molecular system. The interface may have 
a graphical or data file input (or both). The biochem istry components module 52 
translates the modeler input for a particular mathematical model of the molecular system 
and is divided into translation submodules 60, 62 and 64 for mathematical modeling the 
m olecule(s) ? the force fields and the solvent respectively of the system being modeled. 
There are several modeler and bioche mistry components modules available including, for 
example, Tinker (Jay Ponder, TINKER User's Guide . Version 3.8, October 2000, 
Washington University, St. Louis, MO). 

The paragraph beginning on page 1 1, line 32, has been amended as 

follows: 

With the translated physical parameters from the bioche mistry 
components module 52, the physical model module 54 defines the molecular system 
mathematically. At the core of the module 54 is a multibody system submodule 66. The 
analysis module 56, which communicates with the physical model module 54 and the 
visualization module 58, provides solutions to the computational models of the molecular 
systems defined by the physical model module 54. The analysis module 56 consists of a 
set of integrator submodules 68 which integrate the differential equations of the physical 
model module 54. The integrator submodules 68 advance the molecular system through 
time and also provide for static analyses used in determining the minimum energy 
configuration of the molecular system. The analysis module 56 and its integrator 
submodules 68 contain most of the subject matter of the present invention and are 
described in detail below. 

The paragraph beginning on page 12, line 1 1, has been amended as 

follows: 
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The visualization module 58 receives input information from the 
biochem istry components module 52 and the analysis module 56 to provide the user with 
a three-dimensional graphical representation of the molecular system and the solutions 
obtained for the molecular system. Many visualization modules are presently available, 
an example being VMD (A. Dalke, et al. 9 VMD User's Guide , Version 1.5, June 2000, 
Theoretical Biophysics Group, University of Illinois, Urbana, Illinois). 

The paragraph beginning on page 12, line 3 1, has been amended as 

follows: 

The MBS is an abstraction of the atoms and effectively rigid bonds that 
make up the molecular system being modeled and is selected to simplify the actual 
physical system, the molecule in its environment, without losing the features important to 
the problem being addressed by the simulation. With respect to the general system 
architecture illustrated in Fig. 1, the MBS does not include the electrostatic charge or 
other energetic interactions between atoms nor the model of the solvent in which the 
molecules are immersed. The force fields are modeled in the submodule 62 and the 
solvent in the submodule 64 in the bioche mistry components module 52. 

The paragraph beginning on page 13, line 3 1 , has been amended as 

follows: 

The symbol for the vector from one point to another contains the name of 
the two points. Thus, r PQ is the vector from the point P to point Q. A vector 
representing the velocity of a point in a reference frame contains the name of the point 
and the reference frame: N v p . Certain symbols to be introduced later relate two reference 
frames. In this case, the symbol contains the name of two frames. Thus, l C k is the 
direction cosine matrix for the orientation of frame k in frame L This symbol refers to the 
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direction cosine matrix for a typical body in its parent frame. Thus, l C k (j) indicates the 
direction cosine of the actual body j in question. The left and right superscripts do not 
change with the body index. This is also true for the other symbols. An asterisk 
indicates the transpose: H\k) , for example. A tilde over a vector indicates a 3 by 3 
skew-symmetric cross product matrix: vw = v x w . is an i by i identity matrix., and 
Oi. is a zero vector of length i and £. is an i by i zero matrix. 

The paragraph beginning on page 20, line 2, has been amended as follows: 

With the first kinematic calculations described above, the residual 
computation for the Residual Form method can be determined. A detailed description of 
the Residual Form and its application to molecular modeling is found in the previously 
cited co-pending U.S. Patent Appln. No. 10/053354. entitled, "METHOD FOR 
RESIDUAL FORM IN MOLECULAR MODELING," filed of even date. This 



computation fills in two partitions of the vector 



given previously. The first 



,PuJ 

partition is called p q , the kinematic residual, and the second partition is called p u , the 
dynamic residual. The kinematic residual is computed from the difference between a q , 
which is passed-in from the (implicit) integration submodules 66, and the derivative 
computed by each joint: 

q-W(q)u = p q 

The paragraph beginning on page 23, line 24, has been amended as 

follows: 



An important ingredient of this integration process is formation of the 
Jacobian of the differential equations. This is 
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dy 

Since the function /is itself computed by an algorithm rather than by an explicit formula, 
the Jacobian computation represents a substantial amount of work. In the simplest 
approach, the Jacobian can be formed numerically by differencing the derivative routine. 
This is a delicate operation because the quality of the Jacobian is a tradeoff between 
round-off and truncation errors. Typically half the working precision in the result is 
[retained] obtained by choosing a good perturbation size in the difference scheme. In 
practice, though, this is difficult to do. 

The paragraph beginning on page 24, line 10, has been amended as 

follows: 

In general, G, the iteration matrix used in the Newton loop of the implicit 
integrator has the form G-E~aJ, where E is the identity matrix and a is [be] some 
scalar function of the timestep. See the previously referenced U.S. Patent Appln. No. 
10/053.253 . entitled "METHOD FOR LARGE TMESTEPS IN MOLECULAR 
MODELING," filed of even date, for a description of implicit integrators. Changes in 
step size require refactoring G , but not reforming J . Reforming J is needed only when 
the Jacobian is needed at a new state. G is used in a linear subproblem within a Newton 
loop. The following is solved: 

where r{y) is the residual function for that particular implicit integration method. 

The paragraph beginning on page 27, line 16, has been amended as 

follows:: 
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At a high level 4 the residual computation can be considered to depend 
upon two kinds of forces: 'motion forces' and external forces. The motion forces are 
computed directly by the multibody system. The external forces are available to the 
multibody system from a force modeling routine that computes the various interatomic 
forces such as electrostatics and solvents. A similar procedure is followed when 
computing the Jacobian. The multibody system builds the Jacobian of the motion forces, 
and combines it with the Jacobian of the external forces. 

The paragraph beginning on page 27, line 3 1 , has been amended as 

follows: 

The same exchange occurs to compute Jacobians. The native Jacobians in 
their intrinsic coordinates are [be] brought into the MBS coordinates. This requires the 
use of the chain-rule to transform between intrinsic and the MBS generalized coordinates. 
It is important that each effect co-computes its function value and Jacobian, because 
many of the same terms are needed for each computation. Each effect is transformed into 
a set of spatial loads T effect {k) , where k is the index of a generic body in the system. The 

totality of these effects is given the symbol T(k) . 

The paragraph beginning on page 34, line 3, has been amended as follows: 

Note that the row-reduction procedure only needs to be done once before 
computing the Residual Jacobian. The overhead of performing the reduction is more than 
offset by the reduced cost of the smaller matrix vector products which must be formed. 

Note that in forming there is no need to save the elements of the atomic Force 

dr(p,s) 

Jacobian. That is, each element omy nee d s to be available while its contribution 

dr(p,s) 
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to — is being computed. So, more than one element of the big Force Jacobian is 

dr(p 9 s) 

not required at a time. 

The paragraph beginning on page 37, line 17, has been amended as 

follows: 

2. Given (q,u) compute z = -M~ l p u (q 9 u 9 0) using the Direct Form method, [Also] 
Then v set u = z and recompute A(k) , then p u , which recomputes f(k) . 

The paragraph beginning on page 38, line 17, has been amended as 

follows: 

Each interbody direction cosine matrix (and all the joint-specific) quantities 
depend only on the generalized coordinates of an individual joint. Thus, 
l dC k {k) is nonzero only when the derivative is taken with respect to any of the 
coordinates for body L To properly 'seed' the recursions being [studying] 
studied , a vector dq is passed in to the routine. For Jacobian computation we set 

one entry is set to 1, and all the other entries to 0. Then, the needed preliminary 
quantities are generated by a typical loop: 

dq(k) 

The partial derivatives of the direction cosine matrices are generated analytically 
and displayed in the section, "Interface Contraction" above. These partial 
derivatives do not depend upon the particular column of the Jacobian that is being 
computed. Setting a particular entry of dq to 1, and all the rest to 0 has the effect 

of annihilating the correct subset of the seed quantities. 
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9q k 



, the partial derivative of the interbody spanning vector is given by 



dr Q& (k) 

d( lk 
dr Q& {k) 

dr Q ' Qk (k) 
dr 6& (k) 



= A(k), slider 



= 03,4, ball 



= 0,, pin 

^0 0 0 0 1 0 0 

0 0 0 0 0 1 0 

0 0 0 0 0 0 1 



, free 



A(k) here refers to the body's sliding axis which connects it to its parent body. 
dH(k) 



-, the partial derivative of the joint map is 



dH(k) 



= 0%, pin, slider 



dH(k) 

8H(k) 
d <l k = 



d'C*(k) 
3q k . 



for k = 1,2,3,4 free 



= 06 for k =5,6,7, free 



The paragraph beginning on page 39, line 13, has been amended as 

follows: 

Since arbitrary perturbations to a set of Euler parameters do not produce a 
pure rotation, column contraction cannot be used when computing the 
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corresponding column of the Jacobian. The row-reduced Force Jacobian can still 
(and must) beused.[.] 

The paragraph beginning on page 42, line 13, has been amended as 

follows: 



10. Back-solve the [ — ] result of previous step with the mass- 
ed dq 

du 

matrix to obtain the desired — : 

dq 

dq dq 

The back-solve operation is accomplished in the Direct Form method routine by 
processing a residual vector into a u vector. The Second Kinematics Calculations only 
needs to be performed once for the whole Jacobian, since the back-solves are done at the 
nominal value of the state. In fact, the Second Kinematics routine must have been called 
in Step 1 while computing z, so the variables should still be cached. 

The paragraph beginning on page 44, line 7, has been amended as follows: 

14. Back-solve the [— ] result of previous step with the mass- 

du dq 

du 

matrix to obtain the desired — : 

du 

du du 

The back-solve operation is accomplished in the Direct Method routine by processing a 
residual vector into a u vector. The Second Kinematics Calculations only need to be 
performed once, since the back-solves are done at the nominal value of the state. In fact, 
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the Second Kinematics routine must have been called in Step 1 while computing z, so the 
variables should still be cached. 

The paragraph beginning on page 47, line 21, has been amended as 

follows: 

Therefore, while the foregoing is a complete description of some of the 
embodiments of the invention, it should be evident that various modifications, 
alternatives and equivalents may be made and used. Accordingly, the above description 
should not be taken as limiting the scope of the invention which is defined by the metes 
and bounds of the appended claims. 
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